#!/bin/csh -f
##############################################################################
#job file to produce postscript graphic files of tisc results with GMT 4.0
#	Daniel Garcia-Castellanos, 1994-1998.
##############################################################################
#Syntax: 	tisc.gmt.job 'model-root-prj' 
##############################################################################
#Note that:
# -This script uses the bc unix command as a calculator. 
##############################################################################

set prj 	= $1

set width = 8

source $tisc_dir/script/tisc.common.gmt.job

set Region=$Region_plot

#UPPER LEFT: TOPOGRAPHY
psbasemap -JX$size -R$Region -B$ticklabels\:"x (km)":/$ticklabels\:"y (km)":NseW \
	-Y17 -X2 -K -P >! $ps
source $tisc_dir/script/tisc.common.topo+drainage.gmt.job
pstext	-JX -R -O -G0 -W255/255/255 -K <<END >> $ps 
	$xtime	$ytime	13 0 1 9	$Timenow Ma
END


#UPPER RIGHT: GEOLOGY or THIN SHEET VELOCITY FIELD & UNIT THICKNESS
psbasemap -JX$size -R$Region -B$ticklabels\:"x (km)":/$ticklabels\:"y (km)":Nsew \
	-X$horz_shift -K -O >> $ps
if (-r $prj.vel) then 
	#THIN SHEET
	source $tisc_dir/script/tisc.common.thinsheet.gmt.job
else
	#GEOLOGY:
	awk '{if (NR==2) for (i=3; i<=NF; i++) dens[i]=$i; if (NR>2) {for (i=NF; i>3; i--) if ($i-$(i-1)>10) break; if ($NF>=0) print $1,$2,dens[i]; else print $1,$2,1000}}' \
		$prj.hrz | xyz2grd -G$tmp.geol.grd.tmp -I$dx/$dy -R
	source $tisc_dir/script/tisc.common.geol+sedload.gmt.job
endif


#LOWER LEFT: 
psbasemap -JX -R$Region -B$ticklabels\:"x (km)":/$ticklabels\:"y (km)":nseW \
	-X-$horz_shift -Y$vert_shift -O -K >> $ps
if (-r $prj.st) then
	#SEDIMENTS
	source $tisc_dir/script/tisc.common.sediment_thickness.gmt.job
	if (-r $prj.xyzt) then
	    awk '{print $1,$2,$(NF-1);}' $prj.xyzt | \
		xyz2grd -G$tmp.subs.grd.tmp -I$dx/$dy -R$Region -H2
	    grdcontour $tmp.subs.grd.tmp -JX -R -C100  -A1000f7/255/255/255 -G5 -W1/0/0/0 \
		-L0/2000  -O -K -T+0.5c/0.1c:UD >> $ps
	    grdcontour $tmp.subs.grd.tmp -JX -R -C1000 -A1000f7/255/255/255 -G5 -W3/0/0/0 \
		-L0/100000  -O -K -T+0.5c/0.1c:UD >> $ps
	    grdcontour $tmp.subs.grd.tmp -JX -R -C100  -A1000f7/255/255/255 -G5 -W1/0/0/0ta \
		-L-2000/0 -O -K -T-0.5c/0.1c:UD >> $ps
	    grdcontour $tmp.subs.grd.tmp -JX -R -C1000 -A1000f7/255/255/255 -G5 -W3/0/0/0ta \
		-L-100000/0 -O -K -T-0.5c/0.1c:UD >> $ps
	    pstext	-JX -R$Region -O -G0 -K <<END >> $ps 
		#$xtitle	$ytitle	14 0 1 3	Elastic thickness
END
	endif
else
if (-r $prj.xyzt) then
	#DEFLECTION:
	source $tisc_dir/script/tisc.common.deflection.gmt.job
endif
endif



#LOWER RIGHT:
set dx_scale = `echo \-1.5 | bc -l`
if (-r $prj.pfl) then
#2D CROSS SECTION:
	set height_CS = 3.5
	set vert_shift_CS = `echo 0 | bc -l`
	set width_CS = `echo  $width  | bc -l`
	set horz_shift_CS = `echo  $width + $separation  | bc -l`

	psbasemap -JX$width_CS/$height_CS -R0/$dmax/$zmin/$zmax \
		-Ba100f10:"distance (km)":/a2f1:"z (km)":SeW \
		-O -K -X$horz_shift_CS -Y$vert_shift_CS >> $ps

	source $tisc_dir/script/tisc.common.cross_section.gmt.job
else
if (-r $prj.eet) then
#EET:
	source $tisc_dir/script/tisc.common.EET.gmt.job
else
set dx_scale = `echo $width - 1.5 | bc -l`
endif
endif


#PALETAS DE COLORES:
set width_scale = `echo $width \* 2 \- 1.5 | bc -l`
psscale -C$tmp.topo_palet.tmp 	-D$dx_scale/-1.0/$width_scale/.3h -B/:"Elevation (m)": -L -O -K >> $ps
#if (-r $prj.vel) psscale -C$tmp.unit_thick_cpt.tmp 	-D-1.5/-1.25/$width_scale/.35h -B/:"Unit thickness (km)": -O -L -K	>> $ps
#psscale -C$tmp.eet_cpt.tmp 	-D$dx_scale/-4.0/$width_scale/.3h -B/:"EET": -O -L -K	>> $ps
psscale -C$tmp.eros_cpt.tmp 	-D$dx_scale/-2.5/$width_scale/.3h -B/:"Erosion rate (mm/yr)": -L -O -K >> $ps
psscale -C$tmp.sed_cpt.tmp 	-D$dx_scale/-4.0/$width_scale/.3h -B/:"Sediment thickness (m)": -L -O -K >> $ps
#psscale -C$tmp.subs_cpt.tmp 	-D$dx_scale/-7.5/$width_scale/.3h -B/:"Subsidence (m)": -L -O -K >> $ps

pstext -JX -R -O <<END>> $ps
END

rm $tmp.*.tmp 
